home *** CD-ROM | disk | FTP | other *** search
/ Sprite 1984 - 1993 / Sprite 1984 - 1993.iso / src / lib / m / cabs.c < prev    next >
Encoding:
C/C++ Source or Header  |  1988-07-11  |  6.3 KB  |  221 lines
  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms are permitted
  6.  * provided that this notice is preserved and that due credit is given
  7.  * to the University of California at Berkeley. The name of the University
  8.  * may not be used to endorse or promote products derived from this
  9.  * software without specific prior written permission. This software
  10.  * is provided ``as is'' without express or implied warranty.
  11.  *
  12.  * All recipients should regard themselves as participants in an ongoing
  13.  * research project and hence should feel obligated to report their
  14.  * experiences (good or bad) with these elementary function codes, using
  15.  * the sendbug(8) program, to the authors.
  16.  */
  17.  
  18. #ifndef lint
  19. static char sccsid[] = "@(#)cabs.c    5.2 (Berkeley) 4/29/88";
  20. #endif /* not lint */
  21.  
  22. /* HYPOT(X,Y)
  23.  * RETURN THE SQUARE ROOT OF X^2 + Y^2  WHERE Z=X+iY
  24.  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
  25.  * CODED IN C BY K.C. NG, 11/28/84; 
  26.  * REVISED BY K.C. NG, 7/12/85.
  27.  *
  28.  * Required system supported functions :
  29.  *    copysign(x,y)
  30.  *    finite(x)
  31.  *    scalb(x,N)
  32.  *    sqrt(x)
  33.  *
  34.  * Method :
  35.  *    1. replace x by |x| and y by |y|, and swap x and
  36.  *       y if y > x (hence x is never smaller than y).
  37.  *    2. Hypot(x,y) is computed by:
  38.  *       Case I, x/y > 2
  39.  *        
  40.  *                       y
  41.  *        hypot = x + -----------------------------
  42.  *                         2
  43.  *                sqrt ( 1 + [x/y]  )  +  x/y
  44.  *
  45.  *       Case II, x/y <= 2 
  46.  *                                   y
  47.  *        hypot = x + --------------------------------------------------
  48.  *                                       2 
  49.  *                                 [x/y]   -  2
  50.  *               (sqrt(2)+1) + (x-y)/y + -----------------------------
  51.  *                                       2
  52.  *                              sqrt ( 1 + [x/y]  )  + sqrt(2)
  53.  *
  54.  *
  55.  *
  56.  * Special cases:
  57.  *    hypot(x,y) is INF if x or y is +INF or -INF; else
  58.  *    hypot(x,y) is NAN if x or y is NAN.
  59.  *
  60.  * Accuracy:
  61.  *     hypot(x,y) returns the sqrt(x^2+y^2) with error less than 1 ulps (units
  62.  *    in the last place). See Kahan's "Interval Arithmetic Options in the
  63.  *    Proposed IEEE Floating Point Arithmetic Standard", Interval Mathematics
  64.  *      1980, Edited by Karl L.E. Nickel, pp 99-128. (A faster but less accurate
  65.  *    code follows in    comments.) In a test run with 500,000 random arguments
  66.  *    on a VAX, the maximum observed error was .959 ulps.
  67.  *
  68.  * Constants:
  69.  * The hexadecimal values are the intended ones for the following constants.
  70.  * The decimal values may be used, provided that the compiler will convert
  71.  * from decimal to binary accurately enough to produce the hexadecimal values
  72.  * shown.
  73.  */
  74.  
  75. #if defined(vax)||defined(tahoe)    /* VAX D format */
  76. #ifdef vax
  77. #define _0x(A,B)    0x/**/A/**/B
  78. #else    /* vax */
  79. #define _0x(A,B)    0x/**/B/**/A
  80. #endif    /* vax */
  81. /* static double */
  82. /* r2p1hi =  2.4142135623730950345E0     , Hex  2^  2   *  .9A827999FCEF32 */
  83. /* r2p1lo =  1.4349369327986523769E-17   , Hex  2^-55   *  .84597D89B3754B */
  84. /* sqrt2  =  1.4142135623730950622E0     ; Hex  2^  1   *  .B504F333F9DE65 */
  85. static long    r2p1hix[] = { _0x(8279,411a), _0x(ef32,99fc)};
  86. static long    r2p1lox[] = { _0x(597d,2484), _0x(754b,89b3)};
  87. static long     sqrt2x[] = { _0x(04f3,40b5), _0x(de65,33f9)};
  88. #define   r2p1hi    (*(double*)r2p1hix)
  89. #define   r2p1lo    (*(double*)r2p1lox)
  90. #define    sqrt2    (*(double*)sqrt2x)
  91. #else    /* defined(vax)||defined(tahoe)    */
  92. static double
  93. r2p1hi =  2.4142135623730949234E0     , /*Hex  2^1     *  1.3504F333F9DE6 */
  94. r2p1lo =  1.2537167179050217666E-16   , /*Hex  2^-53   *  1.21165F626CDD5 */
  95. sqrt2  =  1.4142135623730951455E0     ; /*Hex  2^  0   *  1.6A09E667F3BCD */
  96. #endif    /* defined(vax)||defined(tahoe)    */
  97.  
  98. double
  99. hypot(x,y)
  100. double x, y;
  101. {
  102.     static double zero=0, one=1, 
  103.               small=1.0E-18;    /* fl(1+small)==1 */
  104.     static ibig=30;    /* fl(1+2**(2*ibig))==1 */
  105.     double copysign(),scalb(),logb(),sqrt(),t,r;
  106.     int finite(), exp;
  107.  
  108.     if(finite(x))
  109.         if(finite(y))
  110.         {    
  111.         x=copysign(x,one);
  112.         y=copysign(y,one);
  113.         if(y > x) 
  114.             { t=x; x=y; y=t; }
  115.         if(x == zero) return(zero);
  116.         if(y == zero) return(x);
  117.         exp= logb(x);
  118.         if(exp-(int)logb(y) > ibig )     
  119.             /* raise inexact flag and return |x| */
  120.            { one+small; return(x); }
  121.  
  122.         /* start computing sqrt(x^2 + y^2) */
  123.         r=x-y;
  124.         if(r>y) {     /* x/y > 2 */
  125.             r=x/y;
  126.             r=r+sqrt(one+r*r); }
  127.         else {        /* 1 <= x/y <= 2 */
  128.             r/=y; t=r*(r+2.0);
  129.             r+=t/(sqrt2+sqrt(2.0+t));
  130.             r+=r2p1lo; r+=r2p1hi; }
  131.  
  132.         r=y/r;
  133.         return(x+r);
  134.  
  135.         }
  136.  
  137.         else if(y==y)          /* y is +-INF */
  138.              return(copysign(y,one));
  139.         else 
  140.              return(y);       /* y is NaN and x is finite */
  141.  
  142.     else if(x==x)            /* x is +-INF */
  143.              return (copysign(x,one));
  144.     else if(finite(y))
  145.              return(x);           /* x is NaN, y is finite */
  146. #if !defined(vax)&&!defined(tahoe)
  147.     else if(y!=y) return(y);  /* x and y is NaN */
  148. #endif    /* !defined(vax)&&!defined(tahoe) */
  149.     else return(copysign(y,one));   /* y is INF */
  150. }
  151.  
  152. /* CABS(Z)
  153.  * RETURN THE ABSOLUTE VALUE OF THE COMPLEX NUMBER  Z = X + iY
  154.  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
  155.  * CODED IN C BY K.C. NG, 11/28/84.
  156.  * REVISED BY K.C. NG, 7/12/85.
  157.  *
  158.  * Required kernel function :
  159.  *    hypot(x,y)
  160.  *
  161.  * Method :
  162.  *    cabs(z) = hypot(x,y) .
  163.  */
  164.  
  165. double
  166. cabs(z)
  167. struct { double x, y;} z;
  168. {
  169.     return hypot(z.x,z.y);
  170. }
  171.  
  172. double
  173. z_abs(z)
  174. struct { double x,y;} *z;
  175. {
  176.     return hypot(z->x,z->y);
  177. }
  178.  
  179. /* A faster but less accurate version of cabs(x,y) */
  180. #if 0
  181. double hypot(x,y)
  182. double x, y;
  183. {
  184.     static double zero=0, one=1;
  185.               small=1.0E-18;    /* fl(1+small)==1 */
  186.     static ibig=30;    /* fl(1+2**(2*ibig))==1 */
  187.     double copysign(),scalb(),logb(),sqrt(),temp;
  188.     int finite(), exp;
  189.  
  190.     if(finite(x))
  191.         if(finite(y))
  192.         {    
  193.         x=copysign(x,one);
  194.         y=copysign(y,one);
  195.         if(y > x) 
  196.             { temp=x; x=y; y=temp; }
  197.         if(x == zero) return(zero);
  198.         if(y == zero) return(x);
  199.         exp= logb(x);
  200.         x=scalb(x,-exp);
  201.         if(exp-(int)logb(y) > ibig ) 
  202.             /* raise inexact flag and return |x| */
  203.            { one+small; return(scalb(x,exp)); }
  204.         else y=scalb(y,-exp);
  205.         return(scalb(sqrt(x*x+y*y),exp));
  206.         }
  207.  
  208.         else if(y==y)          /* y is +-INF */
  209.              return(copysign(y,one));
  210.         else 
  211.              return(y);       /* y is NaN and x is finite */
  212.  
  213.     else if(x==x)            /* x is +-INF */
  214.              return (copysign(x,one));
  215.     else if(finite(y))
  216.              return(x);           /* x is NaN, y is finite */
  217.     else if(y!=y) return(y);      /* x and y is NaN */
  218.     else return(copysign(y,one));   /* y is INF */
  219. }
  220. #endif
  221.